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1 . INTRODUCTION  | 

The  need  for  detailed  terrain  models  in  tactical  simulations  involving  ground  I 

combat  has  resulted  in  considerable  manual  effort  being  applied  to  reading  I 

terrain  heights  from  contour  maps . An  increased  demand  for  terrain  models  has  j 

led  to  the  examination  of  ways  of  making  better  use  of  data  graphics  equipment 
such  as  that  used  by  cartographers (ref. 1)  and  sonic  digitizers  suitable  for  map  ■ 

reading,  i 

There  is  one  obstacle  in  the  way  of  using  the  above  equipment  and  thereby  ! 

significantly  reducing  the  manual  effort  required  to  construct  a terrain  model. 

This  is  the  fact  that  the  programs  in  use  require  terrain  information  at  the  } 

intersections  of  the  grid  lines  of  a rectangular  grid  system,  but  map  reading  j 

equipment  produces  data  in  the  form  of  contours.  This  could  be  dealt  with  by  I 

redesigning  the  tactical  simulation  programs  to  use  contour  information  directly, 
but  it  is  much  more  difficult  to  design  a simulation  which  is  not  based  on  a 
rectangular  grid  and  its  running  time  would  be  considerably  greater.  The 
alternative  is  to  change  the  data,  i.e.  convert  the  contour  information  into 
grid-point  information  by  a suitable  interpolation  procedure  and  leave  the  program 
unchanged.  The  latter  seems  to  be  a more  satisfactory  approach.  There  are 
many  ways  of  producing  a bivariate  interpolation  function  which  is  an  exact  fit 
to  a regular  lattice  of  data  points,  but  an  exact  fit  is  rare  when  the  data  is 
irregularly  spaced.  This  is  the  problem  addressed  in  this  Technical  Report. 

Some  techniques  are  available(ref.2,3)  for  drawing  contours  from  both 
regularly  and  irregularly  spaced  data,  but  these  shed  little  light  on  the  problem 
of  converting  contour  data  to  grid  point  data  which  does  not  seem  to  have  been 
dealt  with  specifically  in  the  literature.  Digitized  contour  data  can,  of 
course,  be  treated  as  irregularly  spaced  data,  an  interpolation  being  carried  out 
at  each  grid  point  using  a weighted  sum  of  a number  of  the  surrounding  data 
points.  These  methods  have  been  discussed  in  reference  4.  They  suffer  from  a 
number  of  disadvantages  for  the  type  of  problem  dealt  with  here.  For  example 
they  do  not  make  use  of  the  intrinsic  regularity  of  the  data.  Furthermore, 
unless  "spot  heights"  (i.e.  degenerate  contours  at  the  apex  of  hills  or  centres 
of  depressions)  are  included  in  the  data  and  are  given  special  weighting  the 
resulting  interpolation  can  flatten  these  areas  by  an  amount  which  is  noticeable. 

This  is  especially  significant  in  "line  of  sight"  studies  where  observers  utilize 
hill  tops.  In  such  situations  small  discrepancies  near  the  top  can  make  large 
differences  in  the  area  of  visibility. 

The  solution  presented  here  attempts  to  use  the  additional  information 
contained  in  the  fact  that  the  data,  although  not  on  a regular  lattice,  is  not 
randomly  positioned  but  organised  in  the  form  of  contours.  It  will  be  shown 
that  a knowledge  of  surface  slopes  will  overcome  the  problem  of  flattening  peaks 
and  troughs . 

The  requirement  can  be  summarised  as  follows:  Find  a bivariate  interpolation 
procedure  which  accepts  and  makes  optimum  use  of  contour  data  to  generate 
terrain  values  on  a regular  lattice  of  points,  so  that  a smooth  continuous 
surface  can  be  fitted  to  all  points. 

The  uniqueness  and  especially  the  minimum  curvature  properties  of  bicubic 
splines (ref .5)  made  them  prime  candi'^ates  for  the  solution  of  this  problem. 

Tbe  treatment  of  bicubic  splines  usually  pre-supposes  data  points  at  the 
intersections  of  a rectangular  mesh.  Since  data  was  not  available  in  this  form 
a new  method  had  to  be  developed.  This  was  done  by  the  procedure  described 
below  resulting  in  an  accurate  and  fast  algorithm  for  producing  a terrain  surface 
satisfying  all  the  requirements  stated  in  the  summary  of  the  problem  given  in 
the  previous  paragraph. 
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2.  THE  DATA  BASE 

Since  contours  are  continuous  and  this  was  to  be  a digital  process  the  first 
step  was  to  decide  on  the  data  base  to  be  used.  This  was  a selection  of  (x,y) 
pairs  on  the  contours,  each  pair  being  assigned  the  value  z of  the  contour  on 
which  it  lay.  The  required  rectangular,  and  regular  grid  system  onto  which  the 
contour  data  was  to  be  transferred  was  superimposed  onto  the  contour  map.  Each 
intersection  of  a contour  and  a grid  line  was  taken  as  a data  point  (see  figure  1) . 
If  in  addition  to  this  information  the  slope  of  the  terrain  was  measured  in  the 
positive  direction  at  the  extremities  of  the  grid  lines,  it  would  be  possible  to 
fit  a unique  linear  cubic  spline  to  the  data  points  along  each  x grid- line  and 
each  y grid- line.  This  defined  a network  of  two  sets  of  splines,  one  set 
running  in  a direction  at  right  angles  to  the  other.  Each  line  from  one  set 
crossed  every  line  from  the  other.  At  the  crossing  point  (x^,  y^)  of  splines 

on  the  grid-lines  x = x^,  Y = Yq  data  values  z(Xg)^z(yg)  of  the  two  splines  are 

not  in  general  identical,  (see  figure  2).  The  first  real  problem  was  therefore 
to  determine  a way  of  choosing  a value  of  z between  z(Xg)  and  zCy^)  which 

preserved  the  continuity  and  smoothness  of  both  curves  while  minimizing  curvature 
over  the  four  new  splines.  This  is  done  in  Section  3. 


3.  A NECESSARY  CONDITION  FOR  THE  INTERSECTION  OF  TWO  SPLINES 
Linear  spline  interpolation  can  be  defined  as  follows (ref .5)  - 


Let  A:  a = xo  < xi  < 


b 


I and  a set  of  real  numbers  {z^l  (i  = 0,1,  M)  be  given.  A function  f(x)  of 

I class  C*  * such  that  f(Xj^)  = z^  and  which  has  given  slopes  po  = f*  (xo)  and 

[ ~ ^ points  is  referred  to  as  a spline  function  and  is 

; denoted  by  S^(f;x).  In  particular  if  the  interpolating  function  f(x)  is  a cubic 

t,  polynomial  the  spline  is  called  a linear  cubic  spline.  Let  L(x;  xo , Xi , ...  x^) 

I denote  the  linear  space  of  all  functions  f(x)  of  class  C^  on  the  interval 

i (^  > X,,)  which  are  equal  to  a cubic  polynomial  in  each  of  the  intervals 

[ (Xi  i*  Xjl  i = 1,2,  . . . M,  i.e.  piecewise  cubic.  Then  the  following  theorem 

' holds (ref. 6)  - 

i Theorem  1.  For  each  set  j z® , z' , ...»  z , Po , p.J  of  values  there  exists 

I exactly  one  f(x)  e L(x;  xo , ...,  x^^)  such  that 

[ f(Xj^)  = z^,  i = 0,1,  . . .,  M,  f'  (xo)  = Po,  f*  (Xj^)  = p^ 

i 

t A continuous  curve  f(x)  has  a radius  of  curvature  p(x)  given  by  the  expression 


P(x) 

I 

[ 


(L 


f (X) 


3/2 


The  reciprocal  of  this  expression  is  the  curvature  at  x. 


When  p is  large. 


* This  implies  that  the  second  derivative  of  the  function  exists  and  is  contin- 
uous . 


- 3 - 
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curvature  can  be  approximated  by  f"  (x) . Since  this  can  be  positive  or  negative 
the  square  of  the  curvature  is  often  used,  as  in  the  following  theorem: 

Theorem  2.  Let  A and  1 be  defined  as  before.  Then  of  all  the  functions 

f(x)  having  a constant  second  derivative  on  the  interval  (a,bl  and 
such  that  f(Xj)  = (i  = 0,1,  ...,  M) ; the  spline  function 

S^(f;x)  with  junction  points  at  the  x^  and  with  Cf;a)  = 

roinimizes  the  integral 


If"  (x)l*  dx 


This  theorem,  due  to  Holladay(ref .7) , is  often  called  the  minimum  curvature 
property. 

Consider  a three  dimensional  Cartesian  coordinate  system  in  which  two  linear 
cubic  splines  are  defined.  One,  on  the  line  y = y^  and  the  other 

S^y(8;y)  on  the  line  x =■  Xq  where. 


Ax: 

a = Xo  < xi  < . . . 

...  <Xm  ' b 

X. 

1 

for  all  i. 

Ay: 

o 

II 

A 

A 

...  <y„  . d 

^i 

for  all  i. 

= po,  sj^^(f;b) 

= 

Pm 

SAy(i;0 

= qo,  S^y(g;d) 

= 

Let  Xg  lie  in  the  segment  (x^  and  y^  be  in  the  segment  (/y.j* 

wish  to  consider  an  additional  joint  in  both  splines  so  that  the  new 


X’ 


®Ayg 


yy) . We 
splines 


a 

c 


xo  < X,  < • • . < x^.l  < Xg  < x^  < . . . < x^ 

yo  < y,  < . . . < y^.i  < yc  < yy  < • • • 


b 

d 


with  the  requirement  that 


= h 


and 

SAx,  = ^Ax  ^^u-1* 


In  particular  note  that 
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(f;  X ,) 

Axg*-  ' u-1 

SAy  (g;  Vl^ 

u 

y A ^ 

Sfiy^ts;  V 


We  now  wish  to  determine  four  functions  - 


VP 

SAytgJ  VP 
S^yCg;  V 


Zi  (x)  on  I X •,  xJ 

U- 1 b 

Z2  (x)  on  [ Xgj  xJ 
Z3  (y)  on  Yq] 

u (y)  on  [ y^  ; y^] 


Satisfying  the  relations 


"•(VP  = ^Ax/^^VP 

b 

zi  (xJ  = SX,^(f;  X^) 

2'<  (y,)  = 

''y) 

"‘^VP  = Vp 

b 

= s^y^(g;  y^.j) 

Z2  (X  ) = Sa  (f  • X ) 

'•  O'*  Axg'^  > U-' 

Z4  Cy„)  = 

Sa  (g;  X ) 

zi  (X(,)  = Z2  (Xg)  = 

23  (Xq)  = 

Z4  (Xq)  = h 

Zi  (Xj,) 

= zi  (Xg) 

Z3  (Xg) 

= zi(yj,) 

and  such  that  the  functional  - 


v(zi(x),  Z2(X),  Z3(y),  Z4(y))  = 


( zt'  (x)]^  dx 


X 

r u , 1 1 
* lZ2 

J V 


(x)]^  dx 


* f [ zV  (y)l^  dy  * f ''  I z'/  (y)p  dy 


is  minimized. 
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The  latter  requirement  is  in  keeping  with  the  minimum  curvature  property  of 
the  previously  defined  sections  of  the  splines  (i.e.  outside  the  segments 

The  problem  is  one  of  determining  an  extremum,  i.e.  a functional  which 
minimizes  the  overall  curvature  in  four  functions,  each  of  which  is  fixed  and 
has  a fixed  slope  at  one  point  and  which  meet  at  the  point  (x^,  y^,  h) 

(h  arbitrary)  with  the  added  provision  of  continuity  of  the  first  derivative  in 
both  the  X and  y directions  in  the  neighbourhood  of  (x^,  y^,  h) . 

A necessary  condition  for  an  extremum  can  be  obtained  from  the  "Calculus  of 
Variations"  (see  reference  8).  A set  of  curves  C = (zi(x),  22 (x),  23 (y),  24 (y)) 
providing  an  extremum  must  satisfy  the  Euler-Poisson  equation  which  for 
functionals  of  the  form  - 

v(z(t))  = F(t,  Z,  z',  ...,  2^"^  dt 

to 


is  the  nth  order  differential  equation  - 


F - 
z 


^ F 
dt  z' 


dx" 


.n 

F r,  + ...  + (-1)"  — F (n)  = 0 

z , n z"-  ' 

dx 


(3) 


where  F^  denotes  the  partial  derivative  of  F with  respect  to  p. 

For  the  functional  (2)  the  set  C must  satisfy  the  four  equations  - 


z. 

1 


"i 


-■j-F^+^F//  = 0 

dx  z . dx  z . 

(i  = 1.2) 

1 1 

-^F/+-^F,/.=  0 

dy  Zj  dy-' 

(i  = 3,4) 

(4) 


obtained  by  repeated  use  of  equation  (3).  The  general  solutions  of  these 
fourth  order  differential  equations  are  the  four  functions  - 

Zi  ^ ^i^**  ^il'  ^i2’  ^i3’  ^i4^  ^ 

z.  = z^(y,  C^2*  S3’  S4^ 

each  of  which  involves  four  arbitrary  constants.  These  will  be  obtained  from 
the  boundary  conditions  and  the  fundamental  necessary  condition  for  an  extremum, 
5v  = 0.  ^ 

Since  the  F are  functions  of  z'.'  alone,  we  have  - 


0 and  F r 
z. 


= 0 (i  = 1,  ...,  4) 


so  that  the  equations  (4)  become  - 


[ 
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^ V H = 0 (i  = 

dx  z^ 

1,2) 

d’  p 

= 0 (i 

= 3,4) 

or 

2f'’^(x)  = 0 (i  = 

1.2) 

(y) 

= 0 (i 

= 3,4) 

The  sclutions  are  therefore  the 

cubic  polynomials 

3 

- 

Zi(x)  = 

j=0 

3 

x^  (i 

= 1.2) 

Zj(y)  = 

J-O 

y^  (i 

= 3,4) 

The  equation  (1)  provide  only  13  independent  relations,  therefore  3 additional 
conditions  are  required  to  solve  for  the  sixteen  coefficients. 

The  four  terminate  or  begin  on  the  line  x = x^,  y = y^,.  By  considering 

the  arbitrary  variations  6z^,  Sz^^  and  the  fundamental  necessary  condition 

for  an  extremum,  namely  6v  = 0,  (see  Appendix  I)  the  following  relations  can  be 
derived  - 


Ji-  p 
■3—  F I » 
dx  zi 


d „ 

d r- 

d c 

- -j—  F . . 

X = Xg-  dx  Z2 

X 

+ -] — F . » 

= X|,+  dy  Z3 

- -r~  F ! 1 

y = Xg-  dy  Z4 

F 

. . 

- F . . 

= 0 

Zl 

X = Xj,-  Z2 

X = X(,+ 

F 

) . 

- F 

= 0 

Z3 

y = Xg-  z-* 

y = Xg^ 

= 0 


r(5) 


Since  F , # 
z. 

1 

written  as  - 


2z^  and  therefore-^  F^» » = Iz^  , the  equations  (5)  can  be 


f 9 t r > 
Z|  (Xp) 


9 9 9 , ^ 9 9 9 , ^ 

z-i  (x„)  + Z3  (yp) 


Zl'  (X^)  - zi'  (Xg)  = 


G-'  '^G-' 

0 


III,  •. 

u (yp) 


= 0 


zV  (Xr)  ■ (yp)  " ® 


(6) 


or  since 


(X) 


2a.2  + 6a.3  x.  z.  (x) 
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z\  (y) 


2a.2  * 6a. 3 y,  z.  (y) 


these  can  be  written  as, 

ai3-a23+a33-a43  = 0 

2ai2  + 6ai3  x^  - 2a2  2 - 6a2  3 = 0 ► (7) 

2a3  2 + 6a3  3 y^  - 2a4  2 - 6043  Xq  ~ 0 

These  provide  the  remaining  3 relations  required  to  uniquely  determine  the  sixteen 
coefficients  These  can  readily  be  shown  (after  arranging  the  a^^^  for 

convenience  in  the  subsequent  reduction)  to  satisfy  - 

[CI4  = ll 

where  - 

4.'  = [ ai  0 ,a2o  ,a.io  ,a4o  ,ai  1 ,a2  1 ,a3  I ,a4 1 ,ai  2 ,a2  2 ,a3  2 ,a3  3 ,ai  3 ,a2  3 ,a4  2 ,a4  3 ] 

b.'  = I zi  (x^_  j) , Z2  fx^) , Z3  (yy_  , Z4  (y^) , z\  (Xy_  j) , Z2  (x^j|) , Z3  (yy_  j) , u (y^)  ,0, ...  ,0] 

and 

r , . .z  .3  n 


''v-1  Vi 


v'  v’ 

V w 


1 

2y  3y^ 

V V 

v -v 

G 

3 2 3 

"g  -V  -V 

"g  -V 

4 

3 2 3 

*G  -V  -V 

V -V 

vS 

’'g 

-v*  -V^ 
^G  V 

1 

1 -1  -1 

1 -1 

■2-g 

3*g-3*g 

-1  -1 

^g 

2 

-2 

«*g-®’'g 

2 

-2  -6Vg 

»(RU-TR-I837(A) 


- 8 - 


This  matrix  was  then  reduced  to  triangular  form  by  elementary  row  operations 
(see  Appendix  II),  the  same  operations  being  applied  to  fe,.  Thus  a set  of 
relations  defining  the  a^^  could  be  written  down.  These  are  shown  in 

Appendix  III,  together  with  a program  for  determining  the  values  of  the  a 


for  any  given  boundary  values. 


ij 


4.  BICUBIC  SPLINE  INTERPOLATION  FROM  CONTOUR  DATA 


The  primary  aim  of  this  Technical  Report  was  to  Establish  the  result  given  in 
Section  3.  There  are  a number  of  ways  of  making  use  of  this  result  to  define  a 
surface  which  is  representative  of  the  surface  described  by  the  original 
contours.  To  make  best  use  of  the  properties  engendered  by  the  process  used  to 
derive  the  lattice  point  data,  piecewise  bicubic  interpolation  should  be  used. 
This  is  described  below. 


4.1  Use  of  the  lattice-point  data  alone 


Consider  a rectangular  area  on  which  a data  base  has  been  established  as 
in  Section  2 and  on  which  the  value  of  z has  been  determined  by  the  method 
of  Section  3,  for  each  of  the  lattice  points  (x^,y^)  i=l,  ...,I-1, 

j=l,  ...,J-1.  And  furthermore  that  there  be  given  values  - 


z(x^,yj)  i = 0,1,  j = 
j = O.J,  i = 

Pij  ■ ‘ j 

qij  - ‘ ‘ 


= 0,  ....  I 

= 0,J 

J 


(See  figure  3 which  displays  the  position  of  the  data  graphically).  Then  it 
has  been  shown(ref.6)  that  given  the  above  values,  the  following  theorem 
holds : 

Theorem  3.  There  exists  exactly  one  piecewise  bicubic  function  z(x,y)  of 
the  form  - 

1+2  J+2 

■ Z E 

m=0  n=0 


where  the  ifi  and  are  piecewise  cubic  and  of  class  C^  on 
m n 

R:  X4)  < X < Xj,  yo  < y < Xj. 

This  theorem  establishes  the  existence  and  uniqueness  of  the  desired  function. 
Furthermore  an  efficient  computational  scheme  for  its  evaluation  has  been 


Within  a given  rectangle  R^^ 


given  in  the  same  paper. 

Xj  J ^ y ^ Xj . the  interpolating  function  2(x,y)  equals  a bicubic  polynomial. 


•J 

I 

m,n*0 


y -^(x  - X.  ,)  (y  - y.  ,) 
'mn  1-1^  j-1 


C j(x,y) 


(9) 
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The  problem  has  been  reduced  to  the  easily  handled  case  of  a regular  mesh 
at  the  cost  of  being  removed  one  step  further  from  the  original  data  since 
the  values  2(x,y)  forming  the  data  base  of  Section  2 are  no  longer  involved 
explicitly. 

4.2  Retaining  the  data  base 

In  order  to  fit  a surface  to  both  the  data  base  points  and  the  lattice 
points  it  would  be  necessary  to  divide  the  X-Y  plane  over  the  region  of  int- 
erest into  rectangular  areas  by  connecting  the  data  base  points 
1 (Xj^.y^),  y^  = y^r  on  grid  line  y = yj^  to  adjoining  grid  lines  / Xj,  1 

y ■ by  means  of  perpendiculars,  and  likewise  the  data  base  points 

{{Xi,y^),  Xj^  = x^j,  on  grid  line  x = with  adjoining  grid  lines  x = x^  j, 

X “ X , . The  values  of  z and  either  p or  q at  the  comers  of  the 

resulting  rectangles  could  be  determined  and  hence  linear  cubic  splines  for 

these  connecting  lines  derived.  Also  the  values  of  z and  p or  q at  the 

intersections  of  such  connecting  lines  could  be  obtained.  Some  means  would 

then  have  to  be  devised  to  obtain  the  s.  . at  the  comers  of  each  of  the 

ij 

rectanplcs  formed  by  this  process.  Then  within  each  such  rectangle  a bicubic 
spliae  of  the  form  given  by  equation  (9)  would  be  used  as  the  interpolating 
function.  The  increased  effort  which  would  be  involved  (if  in  fact  a 
suitable  method  for  determining  the  s^j  could  be  devised)  does  not  seem 

warranted  by  the  expected  increase  in  accuracy. 


5.  DISCUSSION  AND  COMPARISON  WITH  INVERSE  DISTANCE  WEIGHTED  INTERPOLATION 

The  existing  approaches  to  the  problem  of  interpolating  from  irregularly 
spaced  data  have  been  summarized  in  reference  4.  The  most  successful  one  is 
based  on  assigning  a weighted  sum  of  local  values  to  the  point  at  which  inter- 
polation is  required.  The  weighting  function  may  be  a pure  inverse  distance 
function,  where  in  the  neighbourhood  of  a data  point  the  interpolated  value  is 
computed  from  the  difference  of  two  almost  equal  numbers,  which  results  in 
significant  computational  error.  Another  disadvantage  for  terrain  modelling 
work  is  the  fact  that  the  method  imposes  zero  directional  derivatives  at  each 
data  point;  a suitable  choice  of  exponent  would  alleviate  this  problem  to  some 
extent.  It  has  been  discovered  empirically  that  a value  of  2 for  the  exponent 
is  best.  Modifications  to  improve  the  weighting  function  include: 

(a)  adding  a small  constant  to  the  denominator  to  improve  the  performance 
near  data  points, 

(b)  using  an  exponential  form  to  improve  smoothness  while  using  only  a small 
number  of  local  data  points, 

(c)  including  direction  in  the  weighting  function,  and 

(d)  (of  particular  relevance  to  the  current  application)  attempting  to  correct 
the  undesirable  property  of  zero  directional  derivatives  at  the  data 
points  by  the  introduction  of  slope  into  the  weighting  function. 

The  process  described  in  reference  4 for  implementing  method  (d)  was  somewhat 
arbitary  and  considering  the  importance  of  slope  in  "line  of  sight"  work 
especially  at  hill  tops  where  it  is  critical,  it  was  considered  that  a new 
approach  was  required,  and  this  has  been  provided  in  Section  3. 

To  illustrate  the  difference  between  various  methods,  an  example  taken  from 
reference  3 of  an  analytic  surface  representing  "a  steep  hill  rising  from  a 
plain"  defined  by  the  funciton  - 


^ g-kx-s)*  + (y-5)"i 


z 
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was  used.  Contours  at  intervals  of  0.2  were  calculated  with  a data  base 
extracted  using  the  technique  described  in  Section  2.  The  relevant  calculations 
for  the  bicubic  spline  method  are  described  in  Appendix  IV. 2,  and  the  results  are 
illustrated  in  figure  4. 

Two  interpolation  procedures  based  on  inverse  distance  weighting  functions  were 

_2 

also  used.  The  first  was  the  simple  case  with  w(d)  = d given  by(ref.4). 


if  dj  0 for  all  N data  points 


if  d.  = 0 
1 


where  P is  the  point  at  which  interpolation  is  required 
2^  i = 1,  ...»  N are  data  points 

dj^  is  the  distance  between  P and  the  i^*'  data  point 
N is  the  number  of  data  points 

This  method  is  poor  for  the  current  application  since  it  produces  a dip  in 
place  of  a crest  at  the  critical  point  - the  apex  of  the  hill.  This  flattening 
of  hill  tops  and  valley  floors  which  is  characteristic  of  weighted  inverse 
distance  smoothing  is  alleviated  to  some  extent  by  incorporating  two  quantities 

Aj  and  which  represent  the  desired  slopes  in  the  x and  y directions  at  the  i^^ 

data  point  and  are  weighted  averages  of  the  divided  differences  of  z about  that 
point.  They  are  given  by  the  expressions  (see  reference  4)  - 


A.  = 

I 


d'^(z^  - z^)(Xj  - x^)' 


j = l 


diD.,  D.P 


d:^ 

1 


j=l 


'LL  STp^XP 


i=l 


j = l 


where  d(D^,  D^)  is  the  distance  between  the  i^^  and  data  points.  The  inter- 
polating function  then  becomes  - 


Z(P) 


N N 

= (z^  ♦ Az^)  ^ ^ ^ * 


i = l 


i = l 


If  dj  ^0  for  all  u data  points 


If  d.  = 0 
1 
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whore 


lA.(x  - xp  + B.(y 


V 


V . dj 


and 


V = 0.l[max{z^j  - min{z^lj/(maxlA?  + B?lj^ 


This  method  removes  the  dip  produced  by  the  previous  function  and  achieves  a 
value  of  0.901  at  the  apex.  The  curve  fitted  by  the  bicubic  splines  method 
however  achieves  the  value  0.968  giving  a better  approximation  to  the  desired 
curve  with  a value  of  1.0  at  the  apex.  The  data  used  and  the  values  of  the 
interpolation  functions  at  selected  points  are  shown  in  Table  2. 


6.  CONCLUSIONS 

A method  has  been  developed  which  enables  more  efficient  use  to  be  made  of  the 
data  from  automatic  map  reading  instruments.  This  can  now  be  processed  by  the 
means  outlined  in  this  Technical  Report  into  a form  which  can  be  used  directly 
in  existing  terrain  and  vegetation  models. 

The  technique  for  determining  the  minimiom  curvature  solution  for  the  inter- 
section of  two  splines  has  more  general  applicability.  For  example  it  would 
provide  a good  method  of  estimating  the  values  of  missing  data  points  in  a 
regular  mesh  data  set. 
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NOTATION 

i = 1,  2,  . . . , 4,  j = 0,1,  ...,  3.  The  coefficient 
in  the  i^^  of  the  four  cubic  functions. 


the  vector  (ai o ,a2o ,a3o ,a4o ,ei i ,a2 i .as i ,a4 1 ,aj 2 ,a22 ,a32 , 
a3  3»ai3,a2  3.a42,a43) 

the  vector  (zi (x^_ j) ,Z2 (x^) ,23 (yy_ j) , Z4 (y^) ,zi (Xy_ j) ,22 

23 (yy_ j) » Z4 (y^) , 0 , 0 , ...,  0) 


dF{xi  z.  , 2.  z. ) 
111 


i = 1.2 


3F(yi  z^,  2^  z.) 


i = 3,4 


F . , F » » 
Zi'  z. 


similar  definitions  to  those  for  F 


interchange  the  i^*^  rows 


H.(k) 

II.  .(k) 

i.J 


multiply  every  element  in  the  i^j^  row  by  a non-zero  scalar  k 

add  to  the  elements  of  the  i^^  row  k (a  scalar)  times  the 
corresponding  elements  of  the  j row 

z ( x^) 

9z(x^,y^) 


X , ,x 

u-1  u 


z'  (yp 

3z(x^,y.) 


an  element  in  row  i,  column  j as  modified  at  operation 
number  k 

the  term  in  b,  corresponding  to 
3z(x^,y  ) 


respectively  refer  to  the  three  Cartesian  coordinate  axes 


particular  values  of  x and  y 


the  coordinates  of  the  inserted  joint 
the  segment  containing  the  value  x^ 
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^v-l’  segment  containing  the  value 

zCXi./i)  the  value  of  z at  the  point  (x^,  y^^) 

z^  a particular  value  z(Xj),  or  z(x^,  y^)  for  a fixed  j 

z^  a particular  value  or  z(x^,  yj)  for  a fixed  i 

The  four  functions  to  be  optimized  with 
respect  to  curvature 

z^  the  value  k = 1,2 

or  ZjjC/j)  k » 3,4 


(x) 
Zj  (y) 


i = 1.2  1 
i •=  3,4  J 
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APPENDIX  I 

A PARTICULAR  CASE  OF  THE  FUNDAMENTAL  NECESSARY  CONDITION 

FOR  AN  EXTREMUM 

An  increment  Av  of  a functional  of  the  form  v(z(x))  can  be  written  as  - 


Av  = v{z(x)  8z)  - v(z(x)) 


If  it  is  possible  to  resolve  this  into  two  conqionents  so  that  - 

Av  = L(z(x),  6z)  + M(z(x),  5z)  max  |5z| 

where  L(z(x),  6z)  is  a linear  functional  in  8z 
max  |8zl  is  the  maximum  value  of  8z 

and 

M(z(x),  dz)  * 0 whenever  max  |Szl  0 

then  L(z(x),  6z)  is  called  the  variation  of  the  functional  and  is  denoted  by  Sv. 
If  the  variation  of  a functional  exists  and  if  v takes  on  a minimum  or  a maximum 
along  z » zo (x)  then  8v  = 0 along  z = zo (x) (ref.8) . This  is  the  fundamental 
necessary  condition  for  an  extremum  and  is  readily  extended  to  functionals 
involving  several  independent  functions. 

The  functional  under  consideration  in  the  main  body  of  the  text  is 


* [Z3'(y)l* 

JVi 


(I.l) 


It  can  be  seen  that  it  consists  of  the  sum  of  four  functionals  of  the  form  - 

t*ti 

v(z(t))  » r F(t,  z,  z , z")  dt 

J t-to 

For  functionals  of  this  form  with  one  fixed  and  one  moveable  boundary  point,  the 
fundamental  necessary  condition  for  an  extremum  becomes  (see  reference  8)  - 


iF-z' 


F^,  -z 


F » * ♦z 
z 


d u I 

dt  **z*'  t=ti 


F » 1 1 . . 
z t-ti 


fiz,+lF^»»lt, 


Sz't 


(1.2) 
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Now  using  the  linearity  of  the  variation  of  a functional  and  interchanging  the 
limits  of  integration,  the  variation  6v  of  equation  (I.l)  can  be  written  as  - 


Xr.*ix  /•  , 


Av  = ^ ( z'l ' (x)  I ^ dx  - A 


u-1 


I zV  (x)J^  dx  + A I zV  Cy)l*  dy 


'v-1 


yp+5y 

- A r ^ [ zi'  (y)I^  dy 


so  that,  using  equation  (I. 2) 


Av  = \ (-1) 


i+1 


F - z . F » - 
^ "i 


Z . T f I + Z . "3 — F » » 

1 z.  1 dx  Z.  , 

1 1—1 


i=l 


x=x, 


F » - ^ F 
L_  z.  dx  z._j 


x=Xg* 


* "g  " ' ^ x=Xg*  * \ 


/ 


2.  F / - z.  ^ Av 

1 z.  1 z.  1 dy  z. 

1 1 1- 


y=yr* 


I'  » - • • 

^ h h 


y=y* 


6 Zf,  ♦ I F » » I * ^ 

(.  z.  y=yG*  y(^ 


Now  5 x^  and  8 y^  are  zero  but  5 z^,  6 z'^  (the  derivative  of  the  variation  with 


respect  to  x)  and  8z*  (the  derivative  of  the  variation  with  respect  to  y)  are 

all  arbitrary  so  that  their  coefficients  must  be  identically  zero.  Noting  that 

F / i = 1,2, 3, 4,  are  zero  this  can  be  written  as  - 
"i 


iL-  p 
— F 1 1 
dx  Zi 


X=X-- 


^ F »# 
dx  Zj 


F .r 
Z| 


yi  r% 

+ -r*  P ^ ' 
x=Xg+  ay 


y-ycr 


d p 
“T“  F f f 
dy  Z4 


= 0 


- F » f 

X=X,-  Z2 

If 


^=v 


= 0 


(1.3) 

(1.4) 


Indicates  derivative  from  either  left  or  right  as  required 
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APPENDIX  II 

THE  REDUCTION  OF  { Cl  TO  TRIANGULAR  FORM 

TIjc  columns  of  the  matrix  ( C)  and  the  corresponding  terms  of  the  vectors  ^ 
and  li  were  initially  arranged  to  facilitate  the  reduction  described  in  this 
Appendix.  'Hie  notation  used  for  elementary  row  operations  is  as  follows: 

H.  . interchange  the  i'^  and  j rows 

1 » J 

If  (k)  multiply  every  element  in  the  i^^  row  by  a non- zero  scalar  k 

H .(k)  add  to  the  elements  of  the  i^^  row  k (a  scalar)  times  the 
* » J ^1^ 

corresponding  elements  of  the  j row. 

An  element  in  row  i,  column  j,  as  modified  at  operation  number  k (below)  is 
denoted  by  R^^j^  and  the  corresponding  term  in  fe  by  R^  R^^. 

The  following  list  shows  the  operation  followed  by  the  changed  row  and  the 
corresponding  element  of  b - 

(1)  iU.iC-l):  0,0,0,-l,x^-xi  , 0,0, -y^,  - Xi  ^ , 0,0,0,  - xi  ^ , 0,  -Yq  , 

-z. 

(2)  H,o,2(-1):  0,0,0, -1,0,  x^-X2  ,0,  -y^,  0,  x^^  - X2^  , 0,0,0,  x^^  - Xj^, 

-Yg  • *Yq  ; -Z2 

(3)  H,  ,,3(-l):  0,0,0,-l,0,0,yg-y,  , -y^,  0,0, y^^  - Yi  ^ , y^.^  - Yi  ^ ,0,0,-yg^ , 

■^G^ ; -Z3 

(4)  il<»,4(l):  4 X "0",  Xg-xi  , 0,0,  y2 -y^.x^^ -xi  ^ ,0,0,0, x^^ -xi  ’ ,0,y2^ -y^^  , 

Y2'*-Y(^’;  Z4-Z1 

(5)  H,  0,4(1):  5 X "0",  Xg-X2 , 0,  y2 -Yq  0,  Xg^-X2^,  0,0,0,  Xg^-X2*,  Y2^-Yq^. 

Y2  ^ - Yq^  ; Z4  - Z2 

(6)  H,  1,4(1):  6 X "0",  Yq-Yi  , Yz -Yc.0,0,yg^ -Yi  ^ , Yq^ -Yi \0,0,y2 ^ -y^,^  , 

y2^-yg^;  Z4  - Z3 

(7)  ll|,s(-xi):  1,7  X "0",  -xi^,  0,0,0,  -2  xi  \ 0,0,0,;  zi  - xi  z/ 

(8)  H9,5(xi-Xq):  7 X "0",  yz-y^,  (x^-xi  )^  , 0,0,0, x^,^ -xi  ^ - (x^-xi  )^  xi  ^ , 0, 

Y2^  - y^^  yz^  - ■,  z4-zi  - (X(,-xi)z', 

(9)  H,3,5(-1):  5 X "0",  -l,0,0,2(Xg-x,),-2Xg,0.0,3(Xg^-x,^),-3X(,^0,0;  -z', 

(10)  Hio  ,6  (xz-Xj.) : 7 x "0"  ,y2 -yjj,0,  (x^-xz  )^  ,0,0,0, x^^ -xj  ^ - (x^-xa  ) 3X2  * , Yi^-Yq*. 

Yz^-Yq^;  Z4  - Z2  -(Xq-xz)  z'z 

(11)  H2,6(-xz):  0,1.7  X "0",-X2^  ,0,0,0, -2x2^  0,0;  Z2-X2Z2 

(12)  H,3.*(1):  8 X "0",  2(xg-xi),  2(xz -x^)  ,0,0,3(Xj,^ -xi  * ) , 3(X2^ -Xg^  ) ,0,0; 


ZI 
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(13)  H3,7(-y,);  0,0, 1.7  x "0",  -yi*,  -2y,\4  x "0";  Z3-yi  zj 

(14)  Hu,7(yi-yg):  7 x "0",  y2 -yg.O.O.  (y^-y,  )^ . y(,^-yi^-(yQ-yi  )3yi*  ,0,0, 

72^ -Yq^;  Z4-Z3-(yg-yi  )z3 

(15)  H,4.7(-1):  7 X "0",-l,0,0,2(yg-y,).3(yQ*-yi*),0,0.-2yQ,-3yg^  -Z3 

(16)  H4,8(-y2):  0,0,0,1,10  x "0",  -y2^-2y2^;  Z4-y2  zi 


(17)  H,4.8(l):  10  X "0",  2(yg-y,),  3(y(,^ -y,  ^ ) ,0.0,2(y2 -y^) , 3(y2'-yg^); 

I I 

24  - 23 

(18)  Hi , ,8  (yg-y2 ) : 10  x "0",  (yg-yi  )\yg^ -yi  ^ -(yg-yi  )3yi \0,0,-(yQ-y2)* , 

y2  ^ -Yq  - (yi -Yq) 3y2  ^ ; Z4  - Z3 - (y^-yi ) Z3 - (y2 -y^)  u 

(19)  Hio,8  (yQ-y2)  : 9 X "0",  (Xg-X2)*  , 0,0,0, Xg^-X2^-(Xg-X2)  3X2*,  -(yQ-y2)*, 

y2’-yG^-(y2-yG)3y2* ; Z4  - z2-(xg-x2)z'2  - (y2-yG)zi 

(20)  H9,8(yG-y2):  8 X "0",  (x^-xi)*,  0,0,0,  Xg* -xi  * - (Xg-Xi ) 3xi  * ,0,- (yg-y2  )*  , 

y2^-yG^-{y2-yG)3y2* ; z4-zi -(xQ-xi)z'i-  (y2-yQ)  u 

(21)  Hi5(Ji):  8 X "0",  1,  -1,  0,  0,  3Xg,  -3Xg,  0,  0;  0 

(22)  H,6(Js):  10  X "0",  1.  Sy^,.  0,  0,  -1.  -Sy^,;  0 

(23)  H9,is:  (interchange  9 and  15) 

(24)  Hi3  ,9  (2(xi -Xg))  : 9 x "0",  2(x2 -Xi  ) ,0,0,-3(xi -x^,)*  ,3(X2*-Xg*  )+6(Xq-xi  )Xq, 

rv  t f 

0,0;  Z2  - zi 

(25)  His  ,9  (-(Xq-xi  )*)  : 9 x "0",  (x^-xi  )*  ,0,0,2(xi -x^)*  , (x^-xi  )*  3X(,, 

(y2  * -yG^ ) - (y2 -yG) 2y2 .Y2 ^ -Yq  - (yz -yG) 3y2 * ; 

Z4 -Zl +(Xg-Xl  )z'i  -(y2-yQ)  Z4 

(26)  Hio(l./(Xg-X2)*):  9 x "0",  1 ,0,0,0. (x^* -X2* - (Xq-X2 )3x2* )/ (Xg-X2 )* , 

( (y2 ^ -yG^ ) - (yz -yc^ 2y2 ) / (Xg-X2 )* , ( (y2 * -Yq  ) - (yz -Yq) 3y2 * )/ 
(Xq-X2)*  ; (Z4-Z2'(Xg-X2)Z2-(y2-yG)Z4)/(Xg-X2)* 

(27)  Hi3  ,10  (2(xi -X2))  : 12  x"0",  -3(xi -x^)*  ,3x2*-3xg*+6(xg-xi  )Xq-2(X2-xi  ) 

R101426,  -2(x2-x,  )R101526,  -2(x2-Xi  )R101626;  zi-z'i 


-2(x2-Xi)  R10R26 

(28)  H,s,io(-(X(,-xi)*):  12  x "0".2(xi -x^)*  , (x^-Xi  )*  3Xj,- (x^-xi  )*  R101426, 

(yz  * -yG^-  (yz  -yG)2y2  - (x^-xi  )*  Rioi526,y2*  -y^*  “ (yz  -yG)  3y2* 
-(Xg-Xi  )*R101626;  Z4-Z1  -(x^-xi  )zt  - (yz  -y^,)  zi  - (x^-Xi  )*R10R26 

(29)  Hi  1,15:  (interchange  16  and  11) 

(30)  Hi4,n  (-2(yg-yi)):  11  x "0",  -3(yQ-yi)*,  0,0,2(y2 -yi ) , 3y2*-3yg*  ♦ 2(yg-yi) 

* t f 

3yg;  Z4  - Z3 
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(31)  H,»  ,,  (-(/jj-yi)'):  11  X "0-',2(y,-y^,)^0,0,y2'-yg-2y2{yl-yg)  + (yg-y,  )^ 

yi  (yz  -yg)  3y2  ^ + (y^-yi ) ^ 3y^ ; 24  - 23  - (y^-yi ) 23 ' 

-(y2-yg)zi 

(32)  H,4.,2(3(yg-yi)M:  12  x '’0",3(yg-y,  )^-3(y(,-y,  )^2(y2 -y. ) .SCy^^-y^M  + 

6(yg-yi)yg  - 3(yQ-yi)^;  2^-23 

(33)  il,ft.,2(-2y,-yg)M:  12  X "0",  -2(y, -yg)\2(y, -yg)\ (y^-y. )' - (yg-y2  )^ 

2yg^-2y2^+3y2^  yQ+3yiyQ(yi -2y^)*2iyi  -y^)^  ; 
z4-23-(yg-yi  )23  - (y2-yg)zi 

(34)  H,3(-l./3(x,-Xg)M:  12  x "0",  1,  (Sxj^ -Sx^^ +6(Xg-x,  )Xg-2 (Xj -x,  ) 

R101426/-3(xi  -Xg)^  ,-2(x2-Xi  )R101526/-3(x,  -x^,)^  , 

-2(x2-xi  )R101626/-3(xi-X(,)^  ; ( 2^ -2', -2(x2 -xi  )R10R26l 

/-3(xi  -Xg)^ 

(35)  H,4.,3(-3(yg-yi)M:  13  x "0"  ,-3(y(,-y,  )"  (1+R131434)  , 2y2 -2y, -3(y^-y,  )" 

R131534 , 3y2 " - Sy^^  +6 (y^-y,  ) y^- 3 (y^-y, ) ' 

(1+R131634);  24' -23 -SCy^-y, )^  R13R34 

(36)  Ills  .1  j(-2(xi-Xp_)’):  13  x "0",  (x^^-xi  )^  3Xg- (x^-xi  )^  R101426  - 2(xi-X(,)^ 

R1 31434 , (y2  ^ -y^ ) - (y2 -y^) 2y2 - (x^-xi  )^  R101526-2 (xi  -Xg)^ 
R131534,  (yj  -y^,^  ) - (y^  -y^)  3y2 " - (x^-xi )"  R101626 
-2(xt-Xg)^  R131634;  24 -21  - (x^-xi ) 21  - (y2 -y^) 2^- (x^-xi  )^ 

R10R26  - 2(X|-Xp)^  R13R34 

(37)  Hi6.i3(2(y,-y^p’):  13x"0",  2(yi -y^,)'"  (1+R131434) , (y^-yi )' - (yg-y2 )"  + 

2(yi-yg)'  R131534,2yp’-2y2'+3y2'  yg+3y,  yG(yi-2yg)  + 

2(yi-yg)'‘  (1+R131634);  24-23-(y(,-yi  )23-(y2-yg) 

24+2(yt-y„)^  R13R34 

Vi 

(38)  H,4(l./-3(yj,-yi)^  (1+R131434)*):  13  x "0",  I 2y2 -2y, -3(yg-y,  )*  R131534]  F38, 

I 3y2'  -3y(,'+6(y(,-y,  )yg-3(yg-y, )'  (1+R131634)) 
F38;  I 2'4-Z3-3(yg-yi  )*R13R34]  F38 
* Factor  referred  to  as  F38 

(39)  H, 5 ,14 (-R151436) : 14  x "0",  R151536-R151436.R141538,  R151636-R151436. 

R141638;  R15R36-R151436.R14R38 

(40)  H,6,,4(-R161437):  14  x "0",  R161537-R161437.R141538,R161637-R161437. 

R141638;  R16R37-R161437.R14R38 

(41)  Hi, (1./R151539) : 14  x "0",1,  [ R151636-R151436.R141638] /R151539; 

I R15R36-R151436.R14R381  /R151539 

(42)  M,6 .15 (-R161540) : 15  x "0",  R161637-R161437. R141638-R161540. R151641 ; 

R16R37-R161437,R14R38-R161540.R15R41 

(43)  IliA  (1,/I  R161637-R161437,R141638-R161540.R15164ll  ) 
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This  sequence  of  operations  reduces  [ C]  to  the  partioned  fora  - 


III 

1 

8x8 

1 

-|  ID) 

I 01 

1 16x8 

8x8 

1 

Where  1 1]  is  a 


ID]  = 


-XI 


2x1 


The  terms 

L_ 

ti  , tj  . . , , t9  are  given  by  - 

1 

ti 

I -X2^  - (Xg-X2  ) 3X2  * 1 / (Xq-X2 

= Xq+2x2  = R101426 

t2 

= 

I yi  * -Yq^  - (yi  -yc)  2y2 1 / (x^-  X2 ) ^ 

= -(y(;-y2)*/CxQ-X2)*  « R101526 

t3 

s 

1 Yi ^ -Yq  - (yj -Yq) 3/2*1/ (Xq-X2 ) * 

» R101626 

U 

s 

I 3x2*-3Xq*+6(Xq-xi  )Xq-2(x2-Xi  ) 

R101426]/(3(xi-Xq)*)  » R131434 

ts 

= 

-2(x2 -XI  )R101526/ (-3(x,  -x^)* ) 

= R131534 

s 

-2 (X2 -XI  )R101626/ (-3(xi  -x^)* ) 

= R131634 

s 

l2(y2-yi)-3(yQ-y,)'  R1315341/I- 

3(yG-yi)*(l*R131434)l 

- 

1 3{y2*-yQ*)+6(y^,-y,  )yg-3(yQ-yi ) 

* (1+R131634)1  /I  -3(yQ-yi  )* (1+R131434)1 

t9 

a 

(R151636  - R151436.R141638)/R151539 

1 


matrix,  1 0l  is  a 

zero 

matrix  and  I D] 

is  given  by. 

-2x1* 

■ 

-X2* 

-2X2* 

' 

-VI*  -2yi’ 

-y2* 

•2V2* 

\ 

i 

3X1* 

l 

2X2 

3x2* 

2vi  3vi* 

2V2 

3V2* 

i 

-1 

•=^G 

1 

ti 

*2 

ts 

-1 

1 

1 

-1 

-1 

1 

t4 

«S 

t6 

1 

t7 

t* 

1 

t* 

1 

mL 
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The  vector  becomes  - 


Z,  -X,  Z| 


Zj-Xj  Z2 


23 -Xl  23 


24-/2  24 


2l 


f 

22 


> 

23 


2i 
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[ Z4  - Z2  - (Xq-X2  ) Z2  - (72  -/g)  24  1 / (Xg-  X2  )* 

0 

0 

I zi  -z'l  -2(x2  -X,  )R10R26l  /(  -3(xi  -x^)^  1 

( z; -Z3 -3(yg-y, )^  R13R341  /(  -3(yg-y, )* (1+R131434)1 

(R15R36-R151436.R14R38)/R151539 

(R16R37-R16 1437. R14R38-R161540.R15R41)/R16 1642 


A table  showing  the  operations  which  affect  each  row  in  the  order  in  which  they 
were  used  (last  on  the  right)  is  given  below: 


Row 

Operation  numbers 

1 

7 

2 

11 

3 

13 

4 

16 

5-8 

Not  changed 

9 

1,4,8,20,23 

10 

2,5,10,19,26 

11 

3,6,14,18,29 

12 

Not  changed 

13 

9,12,24,27,34 

14 

15,17,30,32,35,38 

15 

21,23,25,28,36,39,41 

16 

22,29,31,33,37,40,42,43 

23  - 
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APPENDIX  III 

THE  COMPUTING  ALGORITHM  AND  PROGRAMME  FOR  THE  INTERSECTION  OF  TWO  SPLINES 

From  the  reduction  of  [ C]  and  JI2,  given  in  Appendix  II  the  relations  defining 
the  sixteen  components  a^^ ^ of  the  vector  a in  the  equation  [ C]  4 = fe. 

Section  3 of  the  main  text,  can  be  written  down.  These  are,  in  order  of 
calculation,  and  using  the  row  element  terminology  introduced  in  Appendix  II. 


343 

[ R16R37-R161437.R14R38-R161540.R15R41] / 
(R161637-R161437.R141638-R161540.R151641] 

342 

ss 

R1SR41  - tg  043 

aj3 

c 

R14R38  - ts  a43  - t7  a42 

S|  3 

s 

R13R34  - t6  a4  3 “ tj  a4  2 ■ t4  a2  3 

®3  3 

= 

a2  3 + a*  3 “ ai  3 

®3  2 

s 

3yg  843  + a42  - 3yQ  833 

*2  2 

s 

R10R26  - ts  043  ■ ^2  842  - ti  O23 

ai2 

s 

3Xq  823  “ 3Xq  013  +822 

where  the  t.  terms 

1 

are 

defined  in  Appendix  II. 

Writing  p^,  i ■ 1,2, 3,4  for  the  variables  Xj ,X2 , yj , yi  respectively 

*il  “ 'i  ■ ‘i2  ■ ‘13  ^ ' *>2'*'* 

•iO  * "l  ■ Pi  S * ^Pi’  “i3  * Pi’  “i2  P"  " ■ 

During  the  reduction  of  [ C]  to  triangular  form  certain  restrictions  were 
introduced  by  the  fact  that  the  scalars  k in  the  elementary  row  operations 
H.  (K),  H. . (k)  were  required  to  be  non-zero  so  that  we  have  - 


(7) 

Xl 

0 

(8) 

Xi 

*G 

(10) 

X2 

*G 

(11) 

X2 

0 

(13) 

yi 

=9fc 

0 

(14) 

Yi 

^G 

(16) 

y2 

0 

(18) 

^G 

Yt 

(27) 

Xi 

X2 

(38) 

R131434  = 

J 
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(39)  R151436  # 0 

(40)  R161437  # 0 

(41)  R151539  * 0 

(42)  R161540  # 0 

and  from  the  calculation  of  843 

R161637  - R161437.R141638-R161540.R151641  =/=  0 

The  basic  computing  routine  ROUTINE  1 incorporates  checks  for  the  above 
conditions. 


ooooonoDooooooonooooorjonDon 
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SUBROUTINE  RT I NEl (XI, X2, Yl, Y2, XG, YC^Zl, ZIP, Z2, Z2n, Z5^ Z5n, Z4, Z4P, 
1 NFLAG) 


THIS  SUBROUTINE  ACCEPTS  BOHNPARY  VALUES  XI,YI,ZI  AMn  A PARTIAL 
PERIVATIVE  OF  Z(Xt,YI)  1-1,.. 4 AMP  FINPS  INTERPOLATION  FUNCTIONS 
TO  PETERMINE  THE  VALUE  OF  Z AT  THE  POINT  (XG,YG). 

ERROR  FLAGS: 

NFLAG  ERROR 

1 Xl-0 

2 Xl-XG 

3 X2-XG 

4 X2-0 

5 Yl-0 

6 Yl-YG 

7 Y2-0 

8 Y2-YG 

9 X1-X2 

10  R131434-1 

11  R151436-0 

12  R161437-0 

13  R151339-0 

14  R161540-0 

14  R161637-R1C1437.R141638-R161540.R151641-0 

DIF  IS  THE  MAXIMUM  DIFFERENCE  BETWEEN  THE  THEORETICALLY  EQUIVALENT 
VALUES  OF  ZG1,ZG2,ZG3  AND  ZG4. 


IMPLICIT  REAL*8(A-H,0-Z) 

C 

C COMPUTE  POWERS  TO  SAVE  REPEATED  CALCHLATION  OF  THESE  VM.HES 

X12-X1**2 

X22»X2**2 

XG2-XG  — 2 

X13«X12*X1 

X23=X22*X2 

XG3-XG2-XG 

Y12=Y1**2 

Y22-Y2--2 

YG2»YG**2 

Y13-Y12-Y1 

Y23-Y22-Y2 

YG3-YG2*YG 


C 

C-  -TERMS  REQUIRED  FOR  LAST  STAGES  OF  TRIANGULAR  REDUCTION 

F26-1./(XG-X2)**2 

R01426-XG+2.-X2 

R01526«(Y22-YG2-(Y2-YG)-2.*Y2)*F26 

R01626-{Y23-YG3-(Y2-YG)*3.*Y22)*F2C 

R0R26»(Z4-Z2-(XG-X2)*Z2D-(Y2-YG)*Z4P)*F26 

F34-l./(-3.*(Xl-XG)**2) 

R31434-(5.*(X22-XG2)+2.-(XG-X1)*3.*XG-2.*(X2-X1)*R01426)*F34 
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I 


inii»5»(=-pnir;2r.*2.  *(X2-xi)*r'^»! 
ir.i()5ii--noir.?r>*?  .*(x2-xi.  )*r3ii 
ir.i!’;»i  = (7.?n-2i'^-'’ . *(X2-xi  )*r:6R2r  )*f  31) 
i!Mir.r.=--( xr-xi  )**7*3.*xr-(xr-xi  )**?*nnTit'’ri- 
I2.*(xi-xn)**3*n-^iii';ii 
i;'.  I'j  3 f = (Y??-Yr 
l-7.*(X7>Xr)**3*n3153I» 

P.!j’  r<3r>  = Y23-Yr'''  - ( ''?-Yr  ) *3  . *v??-  ( yr-X7  )**?*'' m r^r, 

i-2.*(xi-xn)**3*r3in3ii 

Rfjn3r=-z‘'+7  6-(xr-''l)*7.lP-(''2-Yr)*zii^-(Xf'-xi  )**'’*rnP2r. 
l-2.*(Xl-Xr)**3*n3R5ft 
Rrii»37  = '’.*(Yl-Yr)**3*(R31.i(3li  + 7 . ) 
r:m i.37  = ( Yn-Yi  >**2- (vr-v? )**?  + '>. *(yi-Y'')**3*r3i''.3u 
rin  f)3  7 = 2 . *Yn3-2  . *Y'’3  + 3 . *Y2?*'T  + 3 . *vi  *Yn*  ( ''1  - . *vp, ) + o * ( yi  - vr  ) 

1 + 2.  *(Y?.-Yr.)**3*R31fi3li 

Rrir.7  = 7.'i-Z'-(  Yr-V)*Z’,p-(vo_Yr  >*211''  + '’.  *(''■' -Yr)<-*3*r3nTli 
r'r,  = l,/(3.*(YC-Y7  )**2*(- I.-R31I1 ’.())) 

Rli"'  :.32  = (2.*Y2-2  .*','1-:.  .*(Yr,-''l  )**2*P3ir3it 

i:i;ir'28=(  3.  *Y22-3.  *”02  + 6.  *Yr*(  Yr-'n  ) -3  . * ( Yr- Yl)  **  ?*  ( ''31 T 3f;  ♦ 1 . ) ) 
i:M’.3r  = (Zlir-Z3n-3  .*(  Yn-Yl)**2*n3R34  )*r3p 

KM5.3n-^Rr.ll>3f)-R31i|3r.*nitl'j3S 
I : ( . 1 i)  ii  0 - R r,  1 3 7 - f 'v  n 1 4 3 7 * R M 3 3 n 
R3if;4i  = (nr.ir.3r-R:>i4  36*R4ir'3r  )/!:3i43n 
RS  1!4 1 = ( RS R3f,-RI,  14  3R*R4R3R  ) / R3 1 ’.R 

r 

C ZLRO  PIVIDF  riirCKG - 

r;rLAn=! 

iFCXi.i'f'.n.  iroTOiun 
riFi.An=nrLAf>i 
ir(xi-xr,.ro.n.  K'.oToinn 
rjri  Ar=Mri  ap.+i 
I F(  (x2-xr.)  .ro.n . )r:nTni(’n 
r!ri.Ar.=riri  AH+i 
ir(X2.rn.n.  iroio-'n'' 
r!FI.Ar''fJFI  AP  + l 
iF(Yi.rn.o. irninioo 
i.Ti.Ar=r!Fi.  3r+i 
ir(  (Yi-''r).rn.f’.  M'PTnion 
RFi  Ar=riFi.Ar+i 
iF(Y2.rn.n.  ir.OTPm 
NFI.Ar=MFI.Af;  + l 
IF((Y2-Yr).Fn.n.  )''OTniOO 
rjFLAC=MFLAn+l 
IF(R31434  . CO. .imTOlPO 
MFI.Ar.  = NFLAC  + l 
|r(RD143G.ro.O. inOTOlOO 
tlFLAG^'NFLAO  + l 
ir(R(,i  437  .Fa.n.  icnioioo 
r)Fi.An=riFLAn+i 
I F(Rr.  1330.ro.  0.  IGOTOIOO 
rjrLAC=NFLAn+i 
iF(Rr)i340.[;n.M.)(;nTOin() 


n n oo  n n 
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NFLAG-NFLAG+l 

IF(R61637-R61437*R41638-R61540*R51641.FO.O. )GOT0100 
GOTO 10 

100  WRITE(3,19)NFLAG 
RETURN 

19  F0RMAT( 16) 

COMPUTE  COEFFICIENTS  A(I,J) 

10  A43-(R6R37-R61437*R4R38-R61540*R5R41)/ 
1(R61637-R61437*R41638-R61540*R51641) 
A42-R5R41-R51641*A43 
A23«R4R38-R41638*A43-R41538*A42 
A13»R3R34-R31634*A43-R31534*A42-R31434*A23 
A33— A13+A23+A43 
A32«-3 . *YG*A33+A42+3 . *YG*A43 
A22=R0R26-R01626*A43-R01526*A42-R01426*A23 
A12-A22-3.*XG*A13+3.*XG*A23 
All»Zin-2.*Xl*A12-3.*X12*A13 
A21«Z2D-2.*X2*A22-3.*X22*A23 
A31-Z3n-2.*Y1*A32-3.*Y12*A33 
A41-Z4D“2.*Y2*A42-3.*Y22*A43 
A10=Zl-Xl*Zin+2.*X13*Al3+X12*A12 
A20-Z2-X2*Z2D+2,*X23*A23+X22*A22 
A30=Z3-Y1*Z3D+2.*Y13*A33+Y12*A32 

COMPUTE  VALUES  OF  HEIGHT  AT  (XG^YG) - 

ZG1-A10+A11*XG+A12*XG2+A13*XG3 

ZG2»A20*A21*XG+A22*XG2+A23*XG3 

ZG3-A30+A31*YG+A32*YG2+A33*YG3 

ZG4-A40+A41*YG+A42*YG2+A43*YG3 

FIND  MAXIMUM  DIFFERENCE  BETWEEN  THE  ZGI 

ni»nABS(zni-ZG2) 

n2=nABS(ZGl-ZG3) 

n3-nABS(ZGl-ZG4) 

n4»nABS(ZG2-ZG3) 

D5»nARS(ZG2-ZG4) 


n6»nABS(ZG3-ZG4) 

Dl  FoDMAXKDl, 02, 03,04,05, 06) 
C 

C OUTPUT 

WRITE(3,12)A10,A11,A12,A13, 

1 A20,A21,A22,A23, 

1 A30,A31,A32,A33, 

1 A40,A41,A42,A43 


WRITE(3,12)XG,YG 
12  F0RMAT(4F15.4) 

WRITE(3,12)ZG1,ZG2,ZG3,ZG4,0IF 

RETURN 

END 
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APPENDIX  IV 


AN  EXAMPLE 


IV. 1 An  example:  the  intersection  of  two  splines 


Consider  the  two  splines  (f;x),  S^y.  (g;y)  defined  in  Section  3 of  the 


text.  Given  the  values  - 


^Vl’^G* 

^*u-l*^^ 

44 

0 

II 

11. 

10, 

, 0. 

.2) 

^’'u*  ^G* 

Z2(Xy),  z'j(x^^-) 

= (11.5, 

11. 

12. 

1, 

-0.2) 

< 

1 

= (10.3, 

11. 

12, 

.2, 

0.4) 

(yv.XG.Z4 

(yv).z4 (yy-)) 

= (11.4, 

11. 

11. 

.6, 

-0.11) 

= (11. ,11 

.) 

we  wish  to  find  the  functions  Zi (x) , zj (x) , zsCy),  Z4 (y)  satisfying  the 
conditions  (1)  in  Section  3 of  the  text.  The  situation  is  depicted 
graphically  in  figure  2. 

ROUTINE  1 is  used  to  evaluate  the  coefficients  a. . and  hence  the  values 
of  Zg  calculated  from  each  of  the  four  cubic  equations  zi , Zi , zj  and  Z4 . 


Also  the  coefficients  a^^  are  used  to  compute  the  product  fC]  4 and  the 

resultant  vector  fe  is  compared  with  the  expected  vector  fe,.  These  values, 
calculated  using  an  IBM  370/168  computer  in  both  single  and  double 
precision,  are  shown  in  Table  1.  The  derivatives  ziCx)!  _ , zi  (x)!  _ , 

23(y)l„  V • „ ..  are  also  tabulated.  In  general  the  results  show 


that  double  precision  adequately  satisfies  all  the  conditions  (1)  and 
single  precision  gives  a good  approximation  which  may  be  adequate  for  some 
purposes . 

The  functions  z.  i = 1,2, 3, 4 are  plotted  in  figure  4 using  the  double 


precision  results.  The  unbroken  lines  represent  the  sections  of  the  cubics 
in  the  required  regions  I x^  j,  x J , ( y^  j,  y^l  and  the  dotted  curves  are 

included  for  interest  and  show  these  functions  either  side  of  their  required 
ranges . 


IV. 2 An  example:  the  interpolation  at  a local  extremum 

The  function  described  in  reference  3 as  a "steep  hill  rising  from  a 
plain"  given  by  the  expression  - 


-l(x-5)*  + (y-5)^j 

2(x,y)  * e (IV. 1) 

was  taken  as  an  example  of  a local  maximum  to  illustrate  the  performance  of 
the  form  of  interpolation,  recommended  in  this  Technical  Report. 

The  function  was  considered  to  be  defined  by  contours  at  intervals  of 
0.2,  without  "spot  heights",  (i.e.  no  degenerate  contour  at  the  point 
(5.,  5.,  1.)).  The  sixteen  relevant  points  in  the  data  base,  defined  as 
specified  in  Section  2 are  shown  below  • 
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3.751364  4.042769  4.285279  4.527619  5.472381  5.714721  5.97231  6.268636 

5.0  5.0  5.0  5.0  5.0  5.0  5.0  5.0 

0.2  0.4  0.6  0.8  0.8  0.6  0.4  0.2 

These  values  are  used  to  define  the  two  splines  S^(f;Xg),  S^(g;yQ)  from 

which  the  values  zi  j),  zjCXy)»  *3  C/y  j).  can  be  obtained.  In 

addition  to  the  data  in  the  above  table  the  values  p. , i = 0,m,  j = 0,M 

are  needed.  These  can  be  obtained  from  equation  (IV. 1)  noting  that  z(x,y) 
is  synnetrical  about  the  line  x » 5,  y ^ 5. 

iSil  a 

x-3. 731364  ^ 'y*3. 731364 

dz4  I 

x=6. 268636  ' 'y*6. 268636 

The  (for  k » 1,2, 3, 4)  are  determined  by  using  the  following  relationship 
cited  in  reference  6 - 


= 0.507454 

= -0.507454 


h-i  * Pi.i  * ’ 


K A2?-> 


i-l,...M-l 

(IV.2) 


where 


Ax,  “ X.  , - X.  and  Az^  = z^"*^  - zf 

i 1+1  1 k 1 ^ 

Using  the  notation  t^  = 2(Axj^  ♦ the  above  relation  can 

be  written  as  - 
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where  the  0^  are  given  by  the  right  hand  side  of  equation  (IV. 1). 

The  sequence  of  elementary  row  operations  Hi(ti),  (H.  . ,(-u.),  H. (R.)) 

XgX^l  X X X 

i ■ 2,  ....  6,  reduces  the  matrix  to  triangular  form,  so  that  - 


p\  = (01  - Axi  po)/Ri 


(0i  - u.  i = 2,  ....  5 


06  = ( (P6  - Axj  p7  ) - U6  0 S ) /R6 


The  p^  can  therefore  be  obtained  from  the  sequence  - 


P6  * 06 

p,  = 0*.  - u.  , p.  ,/R.  i * 5,  ....  1 


This  enables  the  to  be  determined  noting  that  = p^^  because  of  the 
symmetry  of  the  function, 

zi(Xu_i)  = P3 

23(yv-l^  * q3  = P3  zi(yy)  = q4  = p4 

The  other  required  quantities,  which  can  be  obtained  directly  from  the  data 
base  are  - 

fVr  “ (4.527619,  5.,  0.8) 

(x...  Xp,  Z2  (x,))  = (5.472381,  5.,  0.8) 

U u 

Vl'  ^*527619,  0.8) 


(Xg*  y^.  Z4(yy)) 


(5.,  5.472381,  0.8) 


This  information  was  used  in  ROUTINE  1 (see  Appendix  III).  The  value 
obtained  for  zCS.,  5.)  was  0.9682.  The  cuirve  (IV. 1)  together  with  the 
interpolating  function  is  shown  plotted  in  figure  5.  The  derivatives  of 
both  functions  are  also  shown  in  the  same  figure. 


■ 
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TABLE  1.  DATA  FOR  THE  INTERSECTION  OF  THE  TWO  SPLINES  OF  EXAMPLE  IV. 1 


mean  z 
var  z 


z,  (X-) 
Zj (x+) 

Z3 (y-1) 
Z3  (y-*-) 


Single* 

Double* 

precision 

precision 

2855.74 

2853.87822 

-804.72 

-804.18803 

75.641 

75.59007 

-2.3626 

-2.3609826 

-942.69 

-937.03649 

231.21 

229.69779 

-18.535 

-18.399545 

0.49121 

0.48718773 

-2487.51 

-2482.5837 

697.28 

695.90655 

-64.748 

-64.619135 

2.0012 

1.9971962 

1310.9 

1308.3311 

-338.65 

-337.979 

29.4284 

29.370486 

-0.852638 

-0.85097414 

Di f f erence 


Difference 

(%) 


-2.6 

0.67 

■0.06 

0.002 


11.7400 

11.7371 

11.7314 

11.7454 


11.7385 

0.0058 


Single 

precision 


1.75721 

1.7560 

-0.73422 

-0.73347 


11.741209563968 

11.741209563955 

11.741209563956 
11.741209563958 


11.741209563957 

0 


Difference 


0.00121 

-0.00075 


Double 

precision 


1.7569590 

1.7569586 

-0.7322041 

-0.7322034 


Difference 


0.0000004 


0.0000007 


) 
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TABLE  2.  DATA  AND  FUNCTION  VALUES  FOR  EXAMPLE  IV. 2 
Data  base  for  interpolation  examples 


X 

y 

z 

4.2852 

5.0000 

0.6000 

4.5276 

5.0000 

0.8000 

5.4723 

5.0000 

0.8000 

5.7147 

5.0000 

0.6000 

5.0000 

4.2852 

0.6000 

5.0000 

4.5276 

0.8000 

5.0000 

5.4723 

0.8000 

5.0000 

5.7147 

0.6000 

5.5053 

4.4946 

0.6000 

5.5053 

5.5053 

0.6000 

4.4946 

4.4946 

0.6000 

4.4946 

5.5053 

0.6000 

Simple  inverse  distance  function 


X 

y 

z 

5.0000 

5.0 

0.7067 

5.1500 

5.0 

0.7108 

5.2000 

5.0 

0.7156 

5 . 3000 

5.0 

0.7373 

5.4000 

5.0 

0.7789 

5.6000 

5.0 

0.6856 

Inverse  distance  function  with  gradients 


X 

y 

z 

5.0000 

5.0 

0.9013 

5.1500 

5.0 

0.8898 

5.3000 

5.0 

0.8563 

5.6000 

5.0 

0.7099 

2 given  at  all 
grid  line 
intersections 


Height  ^ ^ L o Slope 
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